#include "PElement.h"
#include <iomanip>
#include <cmath>

PElement::PElement():Element(){
	this->cond=Matrix::zeros(TRIANGULAR,TRIANGULAR);
	this->gradien=Matrix::zeros(TRIANGULAR,TRIANGULAR);
	this->conv=Matrix::zeros(TRIANGULAR,TRIANGULAR);
	this->fGen=Vector::zeros(TRIANGULAR);
	this->fFlux=Vector::zeros(TRIANGULAR);
	this->fConv=Vector::zeros(TRIANGULAR);
	this->force=Vector::zeros(TRIANGULAR);
	
	fbc1=fbc2=cbc1=cbc2=0;
}

PElement::~PElement(){}

PElement::PElement(const int &label):Element(){
	this->cond=Matrix::zeros(TRIANGULAR,TRIANGULAR);
	this->gradien=Matrix::zeros(TRIANGULAR,TRIANGULAR);
	this->conv=Matrix::zeros(TRIANGULAR,TRIANGULAR);
	this->fGen=Vector::zeros(TRIANGULAR);
	this->fFlux=Vector::zeros(TRIANGULAR);
	this->fConv=Vector::zeros(TRIANGULAR);
	this->force=Vector::zeros(TRIANGULAR);
	
	this->label=label;
	fbc1=fbc2=cbc1=cbc2=0;
}

void PElement::constructCond(const double &rho, const double &cp){
	cond=Matrix::diag(TRIANGULAR,TRIANGULAR,2,1);
	double c=(area()*rho*cp/12);
	cond*=c;
}

void PElement::constructGradien(const double &kx,const double &ky,const double &w){
	Vector b=B();
	Vector c=C();
	for(int i=1;i<=TRIANGULAR;i++){
		for(int j=1;j<=TRIANGULAR;j++){
			gradien(i,j)=w*(kx*(b(i)*b(j))+ky*(c(i)*c(j)))/4*area();
		}
	}
}

void PElement::constructConv(const double &h,const double &w){
	if(exposed==true && cbc1!=0 && cbc2!=0){
		int c1=CB1();
		int c2=CB2();
		conv(c1,c1)=2.0;
		conv(c1,c2)=1.0;
		conv(c2,c2)=2.0;
		conv(c2,c1)=1.0;
		double c=h*length(c1,c2)*w/6;
		conv*=c;
	}else{
		conv=conv;
	}
}

void PElement::constructfFlux(const double &q,const double &w)
{
	if(exposed==true && fbc1!=0 && fbc2!=0){
		int f1=FB1();
		int f2=FB2();
		fFlux(f1)=1.0;
		fFlux(f2)=1.0;
		double c=q*w*length(f1,f2)/2;
		fFlux*c;
	}else{
		fFlux=fFlux;
	}
}

void PElement::constructfGen(const double &ugen, const double &w)
{
	fGen=Vector::iso(TRIANGULAR,1.0);
	double c=ugen*area()*w/3;
	fGen*c;
}

void PElement::constructfConv(const double &h, const double &Tf, const double &w)
{
	if(exposed==true && cbc1!=0 && cbc2!=0){
		int c1=CB1();
		int c2=CB2();
		fConv(c1)=1.0;
		fConv(c2)=1.0;
		double c=h*length(c1,c2)*w*Tf/2;
		fConv*c;
	}else{
		fConv=fConv;
	}
}

void PElement::constructForce(const double &q,const double &h,const double &w,const double &gen, const double &Tf){
	constructfGen(gen,w);
	constructfFlux(q,w);
	constructfConv(h,Tf,w);
	force=fGen-fFlux+fConv;
}

